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Variational  Methods  in  Cavitational  Flow 

P.R.  Garabedian 

New  York  University,  Courant  Institute  of  Mathematical  Sciences 
251  Mercer  Street 
New  York,  NY  10012,  USA 

ABSTRACT 

The  study  of  cavitational  flow  is  formulated  as  a free  boundary  problem  for  the  Laplace  equation  in 
three  dimensions.  Constant  pressure  free  streamlines  are  determined  by  a variational  principle  for 
the  virtual  mass.  Steepest  descent  applied  to  minimization  of  the  potential  energy  suggests  a natural 
iteration  scheme  to  calculate  the  shape  of  the  cavity  bounded  by  the  free  streamlines.  Numerical 
methods  enable  one  to  estimate  the  drag  and  the  geometry  of  the  flow.  Another  version  of  the 
variational  principle  plays  an  important  role  in  plasma  physics  and  the  theory  of  magnetic  fusion. 
Novel  stellarator  configurations  for  a thermonuclear  reactor  have  been  designed  by  running  large 
computer  codes  based  on  these  mathematical  ideas. 

1.  Introduction 

We  shall  be  concerned  with  steady,  irrotational  flow  of  an  incompressible  fluid  with  a free  surface  on 
which  the  pressure  is  constant.  This  is  a difficult  mathematical  problem  because  the  shape  of  the 
free  boundary  must  be  calculated  as  part  of  the  solution.  Our  approach  is  to  apply  the  principle 
of  minimum  virtual  mass,  in  which  the  answer  appears  as  the  solution  of  a problem  in  the  calculus 
of  variations  that  is  easier  to  treat  both  theoretically  and  numerically.  The  appearance  of  potential 
energy  in  the  analysis  may  be  unconventional  in  fluid  dynamics,  but  we  shall  show  how  this  can  be 
used  to  find  the  drag,  and  afterwards  we  shall  present  another  application  of  the  same  method  in 
magnetohydrodynamics  where  the  physics  becomes  more  natural. 

2.  The  Riabouchinsky  model 

Steady  irrotational  flow  of  an  incompressible  fluid  in  the  plane  is  governed  by  a complex  potential 
C = (j)  + whose  real  and  imaginary  parts  satisfy  the  Cauchy-Riemann  equations 

4*x  — Ipy  i 4*y  — Ipx  • 

The  horizontal  and  vertical  components  of  the  velocity  are  found  from  the  derivative 

w = u — iv  = dC,/dz 

of  £ as  an  analytic  function  of  the  complex  variable  z = x + iy,  and  the  pressure  p is  found  from 
Bernoulli’s  law 

9 

q " p 

— + — = const.  , 

2 p 

where  q = |w|  is  the  speed.  Thus  along  any  streamline  ip  = 0 bordering  a cavity  we  arrive  at  a free 
boundary  condition  of  the  form 

p = const.  , q = const. 

In  the  Riabouchinsky  model  of  cavitational  flow  the  obstacle  consists  of  two  symmetrically  located 
vertical  plates  joined  by  a pair  of  free  streamlines.  A solution  is  obtained  by  mapping  the  upper  half 
of  the  £-plane  conformally  onto  a semicircle  in  the  hodograph  plane  slit  along  a segment  of  the  u-axis. 
Another  interesting  method  of  finding  the  flow  is  to  study  the  geometry  of  the  analytic  function 

g{z)  = I (’{zf  dz  = z , 


Paper  presented  at  the  RTO  AVT  Lecture  Series  on  “Supercavitating  Flows”,  held  at  the  von  Kdrmdn 
Institute  (VKI)  in  Brussels,  Belgium,  12-16  February  2001,  and  published  in  RTO  EN-010. 
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which  has  boundary  values  indicated  at  the  right  along  a free  streamline  where  q = 1.  This  function 
has  to  do  with  the  forces  on  the  vertical  plates  and  is  suggested  by  the  conservation  of  momentum. 
The  solution  of  the  Riabouchinsky  free  boundary  value  problem  can  be  obtained  by  using  the  Scliwarz- 
Christoffel  formula  to  calculate  the  two  auxiliary  functions 

<h(z)  = z + g(z ) , 'F(z)  = z-g(z)  , 

which  map  the  flow  conformally  onto  polygonal  domains  of  the  complex  plane.  These  classical  methods 
apply  primarily  to  plane  flow  and  play  a significant  role  in  the  study  of  supercavitating  hydrofoils. 
However,  we  do  not  discuss  them  in  detail  here,  but  turn  our  attention  instead  to  three-dimensional 
flows  such  as  occur  in  the  study  of  torpedoes. 


3.  The  principle  of  minimum  virtual  mass 

Let  us  consider  a flow  of  water  past  a body  followed  by  a cavity  fl  in  three-dimmensional  space. 
The  velocity  potential  of  the  flow  is  a harmonic  function  </>  whose  normal  derivative  vanishes  at  the 
boundary,  and  at  infinity  it  has  an  expansion  of  the  form 


with  a coefficient  a related  to  the  virtual  mass.  On  the  free  surface  T bounding  the  cavity  it  satisfies 
the  additional  condition 

(V(/>)2  = const. 

because  the  pressure  is  constant  there.  The  shape  of  the  cavity  has  to  be  adjusted  to  meet  this 
nonlinear  free  boundary  condition,  and  that  is  what  makes  the  problem  hard  mathematically. 


Figure  1:  Geometry  of  the  Riabouchinsky  model  of  plane  flow  around  a finite  cavity  f2  displaying  the 
free  stramlines  T and  two  vertical  fixed  boundaries  C. 

The  kinetic  energy,  or  virtual  mass, 

M = J f J(V<j>-  Vx)2  dV 

is  a Dirichlet  integral  over  the  flow  region.  We  also  introduce  the  volume 

K=  JJJdV- 

which  is  an  integral  over  the  cavity  and  the  body.  An  application  of  the  divergence  theorem  leads  to 
the  remarkable  formula 

47ra  = M + V 

relating  the  virtual  mass  to  the  coefficient  a in  the  expansion  of  the  velocity  potential  at  infinity. 

We  proceed  to  discuss  variational  formulas  for  a and  V under  a suitable  hypothesis  about  smoothness 
at  the  boundary  of  the  flow.  Let  us  shift  each  point  on  the  free  surface  by  an  infinitesimal  distance 
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Sn  along  its  normal  n. 
formula 


It  is  evident  that  the  first  order  perturbation  of  the  volume  is  given  by  the 


SV  = - 


II 


Sn  dS 


where  integration  is  performed  over  the  free  surface. 

Next  observe  that  according  to  the  divergence  theorem,  the  harmonic  function  <f>  minimizes  the  Dirich- 
let  integral  regardless  of  conditions  at  the  boundary  because  its  normal  derivative  vanishes  there.  It 
follows  that  the  perturbation  of  the  coefficient  a under  the  shift  Sn  can  be  computed  with  <p  held  fixed. 
Therefore  we  obtain  an  equally  simply  result 


47 xSa 


Sn  dS 


giving  the  first  variation  of  the  virtual  mass.  This  means  that  a flow  satisfying  the  free  boundary 
condition 

(V<£)2  = A 

can  be  found  from  the  variational  principle 


4t xSa  — \SV  = 0 , 


where  A is  constant.  The  result  is  better  stated  in  the  form 


M — aV  = minimum 


of  an  extremal  problem,  where  a = A — 1 is  the  cavitation  parameter. 

The  theorem  we  have  described,  which  we  refer  to  as  the  principle  of  minimum  virtual  mass,  makes 
sense  only  when  the  cavity  is  finite  so  the  improper  integrals  that  occur  are  convergent.  This  becomes 
the  case  for  the  Riabouchinsky  model  of  cavitational  flow.  In  that  model  two  symmetrically  situated 
plates  C are  connected  by  the  free  surface  F of  a cavity  fh  For  plane  flow  the  geometry  is  illustrated 
in  Fig.  1.  Note  that  the  cavitation  parameter  a is  positive. 

A physical  interpretation  of  the  variational  formula  for  the  virtual  mass  just  states  that  work  equals 
force  times  displacement.  Details  of  a mathematical  proof  are  easier  to  understand  for  the  example 
of  a vortex  ring,  which  will  come  up  again  in  our  discussion  of  magnetohydrodynamics.  Suppose  4>  is 
a harmonic  function  in  a torus  D with  a unit  period  and  a vanishing  normal  derivative,  and  let  O*  be 
the  analogous  function  in  a neighboring  torus  D*  obtained  by  a shift  Sn  of  the  boundary.  If  M now 
stands  for  the  Dirichlet  integral  of  <j)  over  I)  and  M*  stands  for  the  Dirichlet  integral  of  (j>*  over  D*, 
what  we  want  to  show  is  that 

M*  -M  = If  {V<t>)2Sn  dS  , 

where  higher  order  terms  are  neglected  and  the  integration  is  performed  over  the  surface  of  D.  This 
is  true  because  the  Dirichlet  integrals  of  </>*  and  <f>  over  D differ  by  a term  of  second  order  in  Sn.  To 
see  that  we  apply  the  divergence  theorem  to  establish  the  identity 

Iff  -V0)  dV  = I J ^0 f-* ) dS- 1 1 j{<t>*-<t>)  A cjidv  = 0 , 

which  asserts  that  because  of  the  natural  boundary  condition  satisfied  by  4>  its  Dirichlet  integral  over 
D is  stationary.  The  desired  result  now  follows  because 

(V0  + V<t>*  - V<t>f  = (V0)2  + 2 V(j)  ■ (V<?i*  - v<t>)  + - V(j>)2  . 


4.  Calculation  of  the  drag 
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For  special  shapes  of  the  fixed  boundary  the  principle  of  minimum  virtual  mass  can  be  used  to 
calculate  the  drag.  In  the  Riaboucnsky  model  suppose  that  the  free  surface  T is  a sheath  connecting 
a symmetric  pair  of  circular  disks  situated  in  planes  perpendicular  to  the  a;-axis.  Let  us  magnify  the 
whole  configuration  by  the  factor  1 + 77,  where  rj  is  a small  positive  number.  For  this  infinitesimal 
normal  displacement  of  the  boundary  of  the  flow  our  variational  formulas  yield  the  relation 

6M  - a6V  = f / [W)2  - A ]Sn  dS  . 

Along  the  free  surface  the  integrand  vanishes,  so 

SM  - <j5V  = rjh  If  [{Vcj)f  - A]  dS  , 

where  the  integrals  are  now  evaluated  only  over  the  disks,  and  2 h is  the  distance  between  them.  On 
the  other  hand,  by  Bernoulli’s  law  the  integral 

d = \ I j[i+<?-w<f>)2]ds 

over  just  one  disk  is  the  drag.  Since  V and  M have  the  dimensions  of  length  cubed,  and  since  the 
shift  5n  is  defined  by  a magnification,  this  establishes  the  remarkable  relationship 

4 hD  = 3 {aV-M) 


between  the  virtual  mass  and  the  drag  for  Riabouchinsky  flow  past  a circular  disk. 

If  the  formula  we  just  established  for  the  drag  were  generalized  to  the  case  of  two  parallel  plates  of 
elliptical  cross  section,  then  a continuous  transition  could  be  made  from  axially  symmetric  over  to 
plane  flow.  It  is  not  hard  to  see  that  invariance  of  the  drag  under  changes  of  the  eccentricity  of  the 
ellipse  would  then  lead  to  a contradiction.  We  conclude  that  it  is  naive  to  assume  that  the  symmetric 
model  of  Riabouchinsky  flow  always  has  solutions  in  three  dimensions. 

To  compute  the  drag  more  generally  one  needs  to  solve  the  cavitational  flow  problem  numerically.  A 
systematic  scheme  to  accomplish  that  is  suggested  by  applying  the  concept  of  steepest  descent  to  the 
principle  of  minimum  virtual  mass.  This  means  that  to  improve  on  a given  approximation  of  the  free 
surface  one  should  shift  it  along  the  normal  by  an  amount  5n  proportional  to  the  error  (V^)2  — A in 
the  free  boundary  condition.  This  concept  has  been  implemented  in  computer  codes  that  run  quite 
succesfully.  We  shall  discuss  the  matter  in  more  detail  when  we  describe  an  application  to  magnetic 
fusion  in  the  next  section. 

Another  attack  on  the  drag  problem  has  been  developed  using  a stream  function  which  satisfies  the 
partial  differential  equation 


Aif> 


e chp 
V 9y 


0 


with  e = 0 for  plane  flow  and  e = 1 for  axially  symmetric  flow.  For  other  values  of  e the  boundary 
value  problem  still  makes  sense,  and  it  can  be  solved  in  closed  form  in  special  cases.  That  enables  one 
to  estimate  the  drag  for  a circular  disk  by  interpolation  from  known  values.  For  a circular  disk  with  an 
infinite  cavity  the  drag  coefficient  is  found  to  be  Cd  — 0.827.  Similarly,  the  contraction  coefficient  for 
a circular  orifice  in  a plane  wall  is  Cc  = 0.59,  a result  differing  from  the  answer  obtained  in  the  case  of 
plane  flow.  Presently  of  course  many  more  complicated  free  boundary  problems  are  being  investigated 
computationally  by  a variety  of  new  methods. 


5.  Variational  principle  of  magnetohydrodynamics 

Another  enlightening  application  of  the  variational  principle  plays  an  important  role  in  the  plasma 
physics  of  magnetic  fusion.  The  concept  of  a thermonuclear  reactor  is  based  on  fusion  of  hydrogen 
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at  very  high  temperatures  to  form  helium  and  release  neutrons  that  penetrate  a blanket  which  is 
heated  up  to  supply  power.  The  hydrogen  is  ionized  to  become  a plasma  so  hot  that  it  should  not 
encounter  material  walls  and  must  be  confined  by  a magnetic  field.  The  charged  particles  of  the 
plasma  tend  to  follow  lines  of  force  in  the  magnetic  field,  whose  geometry  is  therefore  usually  chosen 
to  be  a torus.  Most  of  the  magnetic  lines  sweep  out  nested  flux  surfaces  rather  than  ergodic  regions 
so  that  the  confinement  of  the  plasma  is  adequate  for  fusion  in  a power  plant.  There  is  an  analogy 
here  with  vortex  rings  and  rotational  flow  of  an  incompressible  fluid  governed  by  the  Euler  equations. 
We  shall  review  the  mathematical  theory  of  equilibrium  and  stability  in  magnetohydrodynamics  that 
is  required  to  address  this  problem. 

We  study  the  equilibrium  and  stability  of  a plasma  by  solving  magnetostatic  equations 

V • B = 0,  J x B = Vp 

analogous  to  the  Euler  equations,  where  B is  the  magnetic  field,  J = V x B is  the  current  density,  and 
p is  the  pressure.  In  the  analogy  the  velocity  of  the  flow  corresponds  to  the  magnetic  field  and  the 
stagnation  enthalpy  from  Bernoulli’s  law  corresponds  to  the  pressure  of  the  plasma.  We  shall  explain 
our  variational  theory  in  the  context  of  the  plasma  problem,  which  will  show  how  islands  whose  flux 
surfaces  are  shaped  like  vortices  can  be  modeled  in  practice  by  current  sheets. 

Let  us  rewrite  the  partial  differential  equations  of  magnetostatics  in  a conservation  form 

V • B = 0,  V • T = 0 


involving  the  Maxwell  stress  tensor  T,  which  is  given  by 

Tjk  = BjBk  — 5jk{B~/2  + p)  , 


where  Sj is  the  Kronecker  delta.  To  avoid  assuming  the  existence  of  partial  derivatives,  we  apply  the 
the  divergence  theorem  and  say  that  B and  p define  a weak  solution  of  the  equations  whenever 


E 


Bk^-dV  = 0 , 

OXk 


.///EE^g 


dV  = 0 


over  any  volume  of  integration,  where  ^1,^2,  V’3  and  are  arbitrary  continuously  differentiable  func- 
tions of  compact  support.  The  simplest  example  of  a weak  solution  is  a magnetic  field  with  just  one 
nontrivial  component  B\  that  is  a nondifferentiable  function  of  £2  and  £3,  while  the  pressure  satisfies 
the  condition  B\/2  + p = const. 

These  formulas  are  comparable  to  the  conservation  of  mass 


■v— ' d)v,k 

^ dxk 


0 


and  the  conservation  of  momemtum 

d dp 

E dxu  uiUk  + pdXj 


0 


for  steady  flow  of  an  incompressible  fluid  without  viscosity  in  hydrodynamics.  Another  application 
of  the  divergence  theorem  establishes  that  any  smooth  surface  of  discontinuity  of  a weak  solution 
of  the  magnetostatic  equations  must  be  a current  sheet  in  the  sense  that  it  is  a flux  surface  of  the 
magnetic  field  across  which  B may  have  jumps,  but  B2 / 2 +p  remains  continuous.  It  is  our  contention 
that  small  magnetic  islands  can  be  modeled  computationally  by  such  current  sheets  when  a finite 
difference  scheme  with  adequate  numerical  viscosity  is  applied.  We  use  the  variational  principle  of 
magnetohydrodynamics  to  show  how  this  concept  has  been  implemented  in  practice.  The  analogy 
with  the  problem  of  cavitational  flow  is  evident. 
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The  variational  principle  of  magnetohydrodynamics  enables  one  to  study  questions  of  toroidal  equi- 
librium and  stability  in  plasma  physics  by  considering  the  extremal  problem 


/ f J[B2/2-p]dV 


= minimum 


for  the  potential  energy  subject  to  appropriate  flux  constraints  on  the  vector  field  B,  which  is  supposed 
to  be  divergence  free.  This  leads  in  a natural  way  to  the  Clebsch  representations 


B = Vs  x V(9 


of  B in  terms  of  flux  functions  s and  6 and  potentials  <j>  and  We  make  a nested  surface  hypothesis 
under  which  s becomes  the  single- valued  toroidal  flux  and  p is  a prescribed  function  of  s.  The  flux 
function  0 is  multiple-valued  on  each  torus  s = const..,  with  periods  related  to  the  rotational  transform 
l = /,(s),  which  measures  how  much  a magnetic  line  turns  in  the  poloidal  direction  during  a complete 
circuit  in  the  toroidal  direction.  For  fully  three-dimensional  equilibria  called  stellarators  no  constraint 
is  imposed  on  i in  the  variational  principle  and  the  potential  </>  has  no  poloidal  period,  whereas  for  an 
axially  symmetric  tokamak  i is  a prescribed  function  of  the  toroidal  flux  s and  (f)  acquires  a nontrivial 
poloidal  period  equal  to  the  net  current. 


Figure  2:  Computational  model  of  fusion  plasma  in  a thermonuclear  reactor  based  on  the  stellarator 
concept  for  magnetic  confinement  of  hot  ions  and  electrons.  Twelve  moderately  twisted  modular  coils, 
half  of  which  are  plotted,  produce  a magnetic  field  whose  strength  has  a desirable  symmetry  displayed 
by  the  color  map  of  the  plasma  surface. 

Straightforward  perturbations  of  the  potential  energy  lead  to  variational  equations  that  are  equivalent 
to  classical  magnetostatics.  Before  a standard  integration  by  parts  they  furnish  a practical  definition, 
essentially  in  conservation  form,  of  what  we  mean  by  a weak  solution  without  derivatives.  For  a 
stellarator  at  zero  (3  = 2 (p)/B2  the  result  reduces  to  Dirichlet’s  principle  for  harmonic  functions,  so 
the  configuration  is  stable.  The  NSTAB  computer  code  has  been  written  to  implement  the  variational 
principle  of  magnetohydrodynamics  numerically  and  perform  an  analysis  of  nonlinear  stability  for 
positive  (3.  It  enables  us  to  capture  small  magnetic  islands  in  a convincing  way  when  we  use  a mesh 
size  for  s larger  than  the  island  width. 

In  stellarators  it  is  convenient  to  renormalize  the  multiple-valued  flux  function  9 and  the  Clebsch 
potential  4>  so  they  become  invariant  poloidal  and  toroidal  angles.  If  one  then  represents  the  magnetic 
field  strength  by  a Fourier  series  of  the  form 

\ 

-^2  = X!  Bmn  COS(TO<?  ~ lH  ~ «?#)  , 
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the  coefficients  Bmn,  which  depend  on  the  toroidal  flux  s,  are  known  collectively  as  the  magnetic 
spectrum.  By  solving  the  magnetostatic  equations  for  J and  taking  a divergence  we  can  show  that 
there  is  a comparable  expansion 


J B 

B2 


= p'E 


niBmn 
n — vm 


cos  (md  — [n 


mi}4>) 


for  the  parallel  current,  which  is  analogous  to  the  swirl  in  fluid  dynamics.  In  three  dimensions  the 
small  denominators  n — cm  are  seen  to  vanish  at  a dense  set  of  rational  surfaces  where  /,  = n/m. 
Therefore  smooth  solutions  of  the  equilibrium  problem  that  have  three-dimensional  asymmetry  do 
not  in  general  exist,  so  that  one  should  only  try  to  construct  weak  solutions.  Correspondingly,  the 
only  smooth  vortex  rings  in  steady  flow  must  be  axially  symmetric. 

The  NSTAB  code  includes  a remarkably  robust  calculation  of  the  coefficients  Bmn.  In  the  spectral 
method  that  has  been  used  it  is  helpful  to  filter  the  Fourier  series  defining  the  solution.  This  is 
especially  true  in  the  case  of  the  parallel  current,  which  requires  evaluation  of  a divergent  series. 
It  turns  out  that  many  aspects  of  equilibrium,  stability  and  transport  just  depend  on  the  magnetic 
spectrum,  together  with  the  profiles  of  pressure  and  rotational  transform.  For  tokamaks  there  is  a 
two-dimensional  symmetry  such  that  only  the  column  Bm o differs  from  zero,  and  similarly  for  straight 
two-dimensional  stellarators  only  the  diagonal  Brnrn  is  present.  Because  of  their  two-dimensional 
magnetic  symmetry,  these  equilibria  have  excellent  transport  properties.  The  problem  of  design  is  to 
exploit  this  theory  to  arrive  at  configurations  for  a fusion  reactor  that  meet  a broad  range  of  physics 
and  engineering  requirements. 


6.  Quasiaxially  symmetric  stellarators 

Fast  computer  codes  with  a three-dimensional  capability  are  an  essential  tool  for  the  design  of  efficient 
stellarators.  Here  we  shall  discuss  in  some  detail  a.  stellarator  called  the  Modular  Helias-like  Heliac  2 
(MHH2)  that  has  quasiaxial  symmetry  characterized  by  very  small  coefficients  Bmn  with  n ^ 0.  This 
was  discovered  recently  by  running  the  NSTAB  code.  It  has  just  two  field  periods  and  the  aspect  ratio 
of  the  plasma  is  only  2.5.  Transport  is  almost  as  good  as  that  in  a tokamak  because  of  the  comparable 
symmetry.  Moderately  twisted  modular  coils  to  generate  the  required  magnetic  field  can  be  wound 
on  a control  surface  surrounding  the  plasma,  whose  shape  brings  the  total  rotational  transform  into 
an  acceptable  range  0.5  > i > 0.4.  Ample  space  is  available  for  the  hardware  requirements  of  a power 
plant  (cf.  Fig.  2). 

We  studied  nonlinear  stability  of  the  MHH2  stellarator  by  looking  for  bifurcated  equilibria  over  the 
full  torus  that  do  not  have  the  helical  symmetry  of  the  solution  with  two  field  periods.  The  most 
dangerous  mode  we  found  running  the  NSTAB  code  has  a complicated  ballooning  structure  that  is 
localized  in  an  outer  rim  of  the  plasma.  Equilibrium  seems  to  impose  a more  severe  limitation  on  [3  for 
the  MHH2  than  stability.  Only  for  broad  pressure  profiles  do  good  magnetic  surfaces  prevail  when  (3  is 
raised  as  high  as  5%.  Thus  calculations  by  the  NSTAB  code  used  to  invent  the  MHH2  establish  that 
its  equilibrium  and  stability  limits  on  average  (3  for  a relatively  broad  pressure  profile  lie  somewhere 
between  4%  and  5%.  The  configuration  is  sufficiently  robust  to  allow  for  further  optimization.  Less 
painstaking  computations  of  linear  or  local  stability  that  ignore  the  difficulties  with  parallel  current, 
at  magnetic  resonances  give  misleading  predictions  about  performance.  Numerical  simulations  built 
around  the  construction  of  weak  solutions  seem  to  be  the  best  that  mathematical  theory  can  contribute 
to  the  problem  of  design  in  this  situation. 
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